Initial data for neutron star binaries with arbitrary spins 
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The starting point of any general relativistic numerical simulation is a solution of the Hamilto- 
nian and momentum constraints that (ideally) represents an astrophysically realistic scenario. We 
present a new method to produce initial data sets for binary neutron stars with arbitrary spins and 
orbital eccentricities. The method only provides approximate solutions to the constraints. However, 
we show that the corresponding constraint violations subside after a couple of orbits, becoming 
comparable to those found in evolutions of standard conformally flat, helically symmetric binary 
initial data. We evolve in time three data sets, corresponding to binaries with spins aligned, zero 
and anti-aligned with the orbital angular momentum. These simulations show the orbital "hang-up" 
effect previously seen in binary black holes. Additionally, all three show orbital eccentricities up to 
one order of magnitude smaller than those found in helically symmetric initial sets evolutions. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 97.60.Jd 



I. INTRODUCTION 

Binary neutron stars (BNS) are currently one of the 
most studied objects in astrophysics due to their poten- 
tial as engines for short gamma-ray bursts [l|, and as 
generators of detectable gravitational waves (GW) 
Since detection rate estimations for the advanced inter- 
ferometric detectors are in the range of ~ 0.4 — 400 BNS 
events per year Q, numerical modeling of the last few 
orbits and merger of such systems are essential for the 
interpretation of the corresponding GW signatures. 

Every numerical simulation has a starting point that 
is, essentially, a snapshot of all the fields (gravitational, 
hydrodynamical, electromagnetic, etc.) at a given time. 
Depending on the characteristics of the modeling formal- 
ism, these fields can either be freely specified or con- 
strained by set of conditions. Numerical simulations 
in general relativity that are based on "3+1" -type for- 
malisms are of the latter kind: the fields have to be so- 
lutions of the Hamiltonian and momentum constraints 
to be consistent with the full set of the Einstein field 
equations 65]. The Hamiltonian and momentum con- 
straints are four coupled second order elliptic PDEs that 
are solved numerically through some iterative procedure 
that starts with an initial guess and loops around the 
equations, correcting the fields until some predetermined 
convergence criteria is reached. Since these four equa- 
tions are not enough to determine the ten independent 
components of the spacetime metric, the modeler has 
the freedom to choose additional constraints/conditions. 
When it comes to finding initial states for BNS in circular 
orbits, the most popular approach is the Wilson-Mathews 
conformal "thin-sandwich" scheme @ which consists 
of restricting the solutions with three extra conditions: 
that the spatial 3-metric jij be conformally flat, that the 
slicing be maximal (tr(Kij) — 0, where Kij is extrinsic 
curvature), and that the spacetime be helically symmet- 
ric (simply put, that the fields be time independent in the 
frame that corotates with the binary). The first two con- 
ditions reduce the number of unknowns from ten to five 



(the conformal factor, the lapse function and the three 
components of the shift vector). The third, however, is 
related to another concern of the initial data (ID): the 
need for it to represent "astrophysically realistic" scenar- 
ios. 

A great deal could be discussed about what constitutes 
an astrophysically realistic BNS; particularly since we are 
largely ignorant of the state of matter inside a neutron 
star. However, there are two aspects of these systems 
that most researchers agree on. One is related to the cir- 
cularity of the orbits of a BNS system that has evolved 
in isolation. In cases like these, it is expected that any 
initial eccentricity the binary may have acquired at birth 
will be minimized by emission of gravitational radiation 
7] (for instance, the BNS known as PSR 1913 + 16 is sup- 
posed to end its life with eccentricities of about 10~ 6 [1]). 
Recent estimates of Advanced LIGO event rates indicate 
that the fraction of BNS with eccentricities larger than 
0.01 is, for the most optimistic scenario, below 2% Q. 
While the helical symmetry condition seeks to approxi- 
mate this by demanding exact circular orbits, it actually 
produces non-negligible eccentricities that, as it is shown 
below, surpasses 0.01 by a factor of several (see also [101]). 

The other aspect of astrophysically realistic ID sets is 
that they should, in principle, be able to describe spin- 
ning stars. Finding ways to construct ID for binaries 
with spinning NS has been a problem more difficult to 
address. Neutron stars obey hydrodynamical equations 
such as the general relativistic versions of the continuity 
and Euler equations and any description of fluids in rota- 
tion should be consistent with them. The original work 
of Wilson and collaborators handled the hydrodynamics 
through a lengthy evolution process that, while not prac- 
tical for the production of ID sets, permitted the imposi- 
tion of arbitrary NS spins by the way of angular momen- 
tum drivers that forced the stars to adopt the desired 
rotations. The impractical nature of this method gave 
rise to a search for simpler techniques. One of the first 
to be considered was the special case when the two stars 
are tidally locked or "corotating" [lj], EH ■ Under coro- 
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tation, the fluid is static in the frame that rotates with 
the BNS and the continuity equation is trivially satis- 
fied. These cases are unlikely to exist in nature since they 
would require fluid viscosities unrealistically high [l3|, [l4| 
but they are still very useful to test new algorithms and 
numerical codes. Corotating solutions were followed by 
solutions with null fluid vorticity ( "irrotational" ) which 
would represent binaries with non-spinning stars. This 
formalism was developed to find solutions for BNS with 
(nearly) zero spin by the way of specifying the fluid veloc- 
ity as the gradient of a potential [L5l - [l8l ] . This potential 
is obtained from an additional elliptic equation derived 
from imposing zero vorticity. Since its introduction, this 
formalism has become the preferred method for BNS ID 
production and several groups have developed codes and 
techniques for its implementation (lH Il9l - [26| . All simu- 
lations of BNS in circular orbits performed to date are 
either based on corotating or irrotational ID sets [27j . 
In recent times, several groups have experimented with 
non-conformally flat techniqu es |28l - t30| . Among these, 
the works by Anderson et al. 13111 . Gold et al. [32|], East 
et al. [33( and Kastaun et al. [34j are of relevance to this 
paper and more about them is said in section [TT1 

However, in general, neutron stars in binaries are ex- 
pected to be spinning. A good example of this is the dou- 
ble pulsar PSR J0737-3039 [35[ that could have one of the 
stars spinning at a rate of about ~ 27ms at the time of the 
merger [36| . Numerical schemes to produce ID for spin- 
ning BNS have been presented by Marronetti and Shapiro 
[37| . Baumgarte and Shapiro j3f| and Tichy [H, [Hj]. All 
of these are based on the Wilson-Mathews helically sym- 
metric, conformally flat approximation. Like in the case 
of irrotational BNS, these methods also rely on advanced 
computationally intensive iterative algorithms. 

We present here a method for producing ID for spin- 
ning BNS that also allows for arbitrary orbital and ra- 
dial velocities. This freedom gives more control over the 
orbital eccentricity than in the case of helically symmet- 
ric methods. Our method does not look for solutions of 
the Hamiltonian and momentum constraints: their satis- 
faction is only asymptotic with binary separation. This 
has the advantage of not requiring the numerical solu- 
tion of elliptic equations, thus greatly simplifying its im- 
plementation. To find out the impact of the resulting 
constraint violation, we produced non-spinning BNS ID 
sets, evolved them in time, and compare the results with 
those of simulations starting with irrotational ID gener- 
ated with the LORENE library (MSI- To facilitate the 
simulations, we implemented of our method as a mod- 
ule ("thorn") of the Einstein Toolkit (ET) [HI, SI]. We 
evolved these BNS for about seven orbits before merger 
and show that, for the grid resolutions and binary sep- 
arations studied here, the constraint violations in both 
simulations become comparable after a couple of orbits 
66] . Additionally, we show that our ID set produces or- 
bits that exhibit eccentricities smaller than those result- 
ing from evolving the helically symmetric ID set. Finally, 
we present the evolution of two ID sets with spinning NS 



(spins aligned and counter-aligned to the orbital angu- 
lar momentum) that showcase the ability of our method 
to handle rotating stars. These simulations present the 
orbit "hang-up" effect [43[ in BNS that is also seen in 
0. 

This article is organized as follows. Sections |TT] and 
IIIII introduce our method and present some tests respec- 
tively. In section IIVI we use our approach to construct 
ID for non-spinning and spinning binaries and show the 
results of their evolution in full general relativistic hy- 
drodynamics using the Einstein Toolkit [4l[. Finally, we 
will briefly summarize our findings in SecfVl 

II. CONSTRUCTION OF INITIAL DATA SETS 

Single rotating neutron stars 

Our approximation to ID for BNS starts with the solu- 
tion of an isolated rotating NS in equilibrium. The study 
of stationary rotating NS has been undertaken by a num- 
ber of groups in the past (see [44j and there in). These 
studies have lead to the creation of publicly available 
codes specially designed to find numerical solutions in 
the framework of the theory of general relativity. One of 
such codes is RNS, developed by Stergioulas and Fried- 
man [45|, Rt| and based on the Komatsu-Eriguchi-Hachisu 
method [471 148| which describes the geometry of station- 
ary and axisymmetric rotating NS with a metric of the 
form 

ds 2 = - e A+B dt 2 +e 2C (dr 2 +r 2 d6 2 ) 

+e A ~ B r 2 sin 2 9 (d(f> — D dt) 2 , (1) 

where the metric functions A, B, C and D depend on r 
and 6. The equations for the gravitational and matter 
fields are then solved using a combination of integral and 
finite-differencing techniques |48j |. 

Our initial data 

Our ID is produced following these three steps: 

(i) Calculate the fields corresponding to two isolated 
rotating neutron stars using RNS, 

(ii) Rotate (if needed) and boost independently each 
solution and map them into a single inertial frame, 
making sure that the BNS total linear momentum is 
zero, and 

(iii) Superpose the fields as indicated below. 

Step (i) is straightforward and generates two station- 
ary solutions for rotating NS. Each one is originally given 
in the reference frame of the RNS code: x ,fi . RNS pro- 
vides solutions in polar coordinates with the NS rotating 
around the z' axis. Since we are interested in NS with 
spins arbitrarily aligned, a rotation of the solutions may 
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be required in combination with the boost. To simplify 
the notation, we will ignore here the rotation (i.e., we 
will consider BNS with the spins in the direction of the 
orbital angular momentum). 

Step (ii) is inspired in standard BBH superposition 
methods [ig], [5(| and starts with the Lorentz boost of 
the RNS coordinates into the inertial frame coordinates 



(2) 



where the Lorentz transformation a A^ is a function of 
the corresponding boost velocity v a , with the index a 
labeling each star (a = 1,2). We now map the metric 
functions to the new coordinate system using ([2]) 

a A'(x'n a A!{x v ) 

a B'(x'n a B\x v ) 

a C'(x'n a C'(x v ) 

a D'(x'n a D'(x v ). 

This allows us to write the spacetime metric g'f x as a 
function of the new coordinates x v and Lorentz transform 
it to the inertial frame 



(3) 



where there is no summation over a. Since we are in- 
terested in solutions in terms of a "3+1" decomposition 
of the metric, we extract from a 9T\{x v ) the correspond- 
ing lapse function a a, shift vector a /3 l , and spatial metric 
ajij- Similarly, a mapping/transformation is applied to 
the rest mass and fluid spatial velocity a p' and a v n , to 
obtain the fields a p and a v l . The latter is defined as 
v l = (u l /u° + (3 l )/a, where is the 4-velocity of the 
fluid. 

Step (iii) constructs a single global set of gravitational 
fields by a superposition of these two solutions 



a 



1 



(4) 



Since there is no overlap between the stars, the superpo- 
sition of the hydrodynamics fields is simply 



P = 1P + 2P 
v l = !«' + 2 v l . 



(5) 



In addition to the stellar matter, a pervasive atmosphere 
is added outside the stars with a relatively low density 
Patm = W~ 7 Pmax(t = 0) and zero velocity. The remain- 
ing hydrodynamical fields (pressure and internal energy) 
are set by the equation of state (EOS). Finally, the ex- 
trinsic curvature is calculated using 



Ki 



1 

~2a 



(6) 



where £p is the Lie derivative along the direction of the 
shift vector. 

Superpositions such as this have been used for BNS 
simulations in the past. In particular, Anderson et 



al. [31(, Gold et al. 1321 a nd, more recently, East et al. 33 1 
and Kastaun et al. [34| have generated ID and evolved 
BNS using either a superposition similar to the one de- 
scribed above or even adding the extra step of act ually 
solving the Hamiltonian and momentum constraints [67| . 
However, ID sets constructed in this way present two un- 
desirable features which are evident during their time 
evolution. The first consists of large oscillations of the 
stellar shape that can be observed by monitoring the 
central density. Figure [1] shows the evolution of the nor- 
malized central density for three different simulations. 
The solid line corresponds to the evolution of a refer- 
ence BNS ID with initial separation of 60 km generated 
by the LORENE libraries, details of which are given in 
Table fl] The dashed line corresponds to a comparable 
ID set (BNS60nl) generated using the steps mentioned 
above. BNS60nl presents large oscillations that contrast 
with the flatness of the evolution of LORENE60 (this 
set is also known as G2_I12vsl2_D5R33_60km). Details of 
these simulations are given in the next section. While 
the amplitude of the oscillations is partially due to the 
grid sparseness (as detailed in the next section, the nu- 
merical grids are at the low end of the resolution range 
currently used in the field), the difference between both 
runs indicates that there are additional causes at play. 
One of them is the fact that simple superposition does 
not make any attempt at coupling the hydrodynamical 
fields (obtained from single NS solutions) with the met- 
ric fields given by These oscillations will diminish 
with increasing binary separation (see Table I in [331 ]) 
and could, to some extent, be controlled by introducing 
fluid viscosity terms in the Euler equations. However, we 
explored alternative ways of minimizing this effect that 
could be easily implemented in the ID generation code. 
Visual inspection of the stellar cross-section on the or- 
bital plane showed that the stars corresponding to the 
LORENE ID set are approximately oval in shape, with 
the diameter along the direction between stellar centers 
larger than the one in the direction of the orbital veloc- 
ity. While our ID also presents this feature (a result of 
the Lorenz contraction in the boost direction), it does so 
at a lesser degree. As an experiment, we tried increasing 
the deformation of the stars in our ID by replacing the 
coordinate transformation ([2]) with 



x v = upv 



(7) 



where the exponent indicates that transformation matrix 
was applied n times. The result of this experiment cor- 
responding to n — 4 is shown as the evolution of the ID 
set BNS60n4 in Fig. Q] We do not offer any justification 
for this "over-contraction" beyond its empirical success. 
Note that this multiple application of the Lorentz ma- 
trix is only performed in the coordinate transformation 
([7]) and not in the transformation of the components of 
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FIG. 1: Maximum rest mass density for a simulation starting 
from a reference ID set (LORENE60 in Table IJ), a superpo- 
sition ID set (BNS60nl) and an "over-contracted" superposi- 
tion ID set (BNS60n4)- M represents the total mass of the 
LORENE60 binary. 



the metric © or the hydrodynamical fields. 

The second weakness of the ID recipe as given is re- 
lated to the orbital eccentricity: using the superposition 
outlined in Eqs. (IHO leads to binaries that exhibit a pro- 
nounced initial eccentricity. Figure [2] shows this effect on 
the behavior of the coordinate separation between stel- 
lar centers for the evolution of the LORENE60 set and 
that of BNS60nl. This problem also diminishes with in- 
creasing binary separation. Again we experimented with 
alternatives to the superposition method that would re- 
duce this spurious eccentricity. We found that rescaling 
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FIG. 2: Coordinate separation for the reference ID set 
LORENE60, the superposition ID set BNS60nl and the su- 
perposition ID set with rescaling BNS60sc. M represents the 
total mass of the LORENE60 binary. 

the hydrodynamical fields using the lapse function has 



the desired effect. 
© with 



We implement this by replacing Eqs. 



a 
a 



Pi 



f>2 



(8) 



This rescaling could be interpreted as a weighted average 
of the hydrodynamics fields by a measure of the space- 
time curvature. The drastic eccentricity reduction caused 
by this rescaling is shown in Fig. [5] as the evolution of 
the set BNS60sc. 

To summarize, our recipe for BNS ID consists in fol- 
lowing steps (i - iii), employing the coordinate transfor- 
mation given in Eq. ([7|) and the calculation of the fields 
given in Eqs. flU), ©, and ©. 



III. NUMERICAL TESTS 

The algorithms described in the previous section have 
been implemented in a numerical module ("thorn") for 
the Einstein Toolkit framework [4l], [42|. The code ac- 
cepts the selection of one or two NS with arbitrary boost 
velocities, spin directions and initial positions. The char- 
acteristics of the NS (mass, spin, EOS) are selected 
through free parameters in RNS which provides the stel- 
lar profiles to be mapped into the ET grid. To test 
our code, we generated ID for single and binary NS and 
evolved them using ET. All the simulations have the same 
grid domain: 5 levels of mesh refinement, provided by 
Carpet thorn [5l|, with box sizes 320M Q (outer bound- 
ary), 120M o , 60M q , 30M q and 15M Q , resulting in a 
resolution of approximately 36 points across the stellar 
diameter. Since it is not our intention to produce high 
quality BNS models but to test the viability of our ID, 
we opted for a computationally affordable and expedient 
numerical setup. The BNS simulations for this article 
make use of xy-plane reflection and 7r-rotation symme- 
tries, since they pertain to systems with identical NS 
with non-precessing orbits. We employ a polytropic EOS 
(p = Kp T ) with K = 123.6 and T = 2.0 for both, the ID 
sets and the evolution. Spacetime evolution is obtained 
through the BSSNOK formalism [Hj], provided by 
the McLachlan/MLJBSSN thorn [EJ]. The general relativis- 
tic hydrodynamic equations are evolved by the GRHydro 
thorn (5f| . We use a Marquina Riemann solver with PPM 
reconstruction of the primitive variables. Kreiss-Oliger 
dissipation is added to the right-hand-sides of the BSS- 
NOK evolution quantities with the dissipation strength 
parameter tdiss — 0.01. 



Isolated Stars 

We start our tests with four single NS cases: with and 
without spin, and stationary and boosted. While these 
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Jadm[M'£\ 


a 
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t m [ms] 
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LORENE60 


60 


13.45 


1.625 


3.005 


9.716 





0.1238 





n/a 
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BNS60nl 


60 


13.45 


1.744 


3.246 


9.155 





0.1460 





n/a 
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BNS60n4 


60 


13.45 


1.689 


3.144 


8.861 





0.1460 





n/a 
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BNS60sc 


60 


13.45 


1.656 


3.079 


9.613 





0.1460 





n/a 


1 


LORENE80 


80 


18.02 


1.625 


3.011 


10.825 





0.1098 





31.76 




BNS80n 


80 


18.02 


1.625 


3.020 


10.497 





0.1265 


-0.005 


31.20 


4 


BNS80u 


80 


18.02 


1.626 


3.034 


12.142 


0.327 


0.1250 


-0.005 


41.71 


4 


BNS80d 


80 


18.02 


1.625 


3.032 


8.984 


-0.327 


0.1265 


-0.005 


25.15 


4 



TABLE I: Parameters of the ID sets used in the binary simulations. Over-contraction Q and rescaling ([8]) were applied to 
all BNS80 sets. The coordinate separation is given by d, while Mo, Madm and Jadm are the rest mass and ADM mass and 
angular momentum respectively, a is the dimensionless spin parameter for each NS, and and v r are the tangential and radial 
components of the boost velocity. Finally, t m is the time to merger (defined by the peak of the rest mass density maximum), 
M* is the ADM mass of the corresponding LORENE set and n is the exponent of the over-contraction (JTJl . 
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FIG. 3: Normalized maximum rest mass density for four single 
NS test cases that include stationary, spinning and boosted 
stars. 



FIG. 4: Hamiltonian and momentum constraint violation (L2 
norm) for the runs of Fig. [3] 



seem rather trivial since they only involve mapped solu- 
tions of the RNS code into the ET grid, they fulfill two 
purposes: to test the mapping/boost algorithms and to 
provide an apt reference in terms of the constraint viola- 
tions that will be useful when analyzing the simulations 
of section IIVI 

The NS rest masses are Mq — 1.625M Q (compaction 
ratio = 0.14). The dimensionless spin parameter a 
(J/M 2 ) is either zero or 0.33 and the boost velocity v 
is either zero or 0.2c. Figure |3] shows the evolution of the 
maximum rest mass density. We observe a stable evo- 
lution for all four cases, with central density oscillations 
consistent with our grid resolution. Figure [4] shows the 
Hamiltonian and momentum constraint violations for the 
four cases under consideration. 



Binaries 

Figure [5] presents a comparison of several met- 
ric fields between a LORENE generated irrota- 



tional ID (LORENE80 in Table H also known as 
G2_I12vsl2_D5R33_80km) and a non-spinning BNS ID set 
with coordinate separation of 80 km produced with our 
method (BNS 8 On in Table 0}. These plots show that even 
the simple superposition of Eq. (j4]) gives good agree- 
ment. We experimented with higher order alternatives 
that, while reducing the separation between curves, did 
not improve noticeably the results of the next section. 
The agreement of the fluid velocity can also be improved 
by adding small positive spins to the NS. We decided to 
neglect that higher order correction in this paper. Figure 
[S]shows the point-by-point convergence of the BNSSOn at 
t = along the line connecting the stars. For the conver- 
gence test we set the initial data on an extended unigrid 
for the resolutions Ax = (0.5/M©, 0.25/M©, 0.125/M©) 
which correspond to (0.738fcm, 0.369A:to, 0.184fcm). All 
our simulations used a central resolution equal to the low- 
est of the ones shown in Fig. [S] (Ax = (0.5/M Q ). This 
figure shows curves consistent with the fourth-order of 
the finite-difference stencils used for spatial derivatives. 
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FIG. 5: Comparison of several fields from the reference ID set 
LORENE80 and our non-spinning ID set BNS80n. The fields 
are plotted along the line connecting the two NS centers. 



IV. BINARY EVOLUTIONS 

A. Comparison with known initial data 

One of the most important tests of any ID set is pro- 
vided by its evolution. Qualitatively good agreement 
between the fields such as the one presented in Fig. [5] 
could lead to large behavior differences after several or- 
bits. In this section we study the evolution of the non- 
spinning BNS ID set BNS80n and compare it with that 
of LORENE80. All the BNS80 sets in this section were 
constructed using over-contraction ([7]) with n = 4 and 
rescaling (JSJ) . Figure [7] shows three snapshots of the rest 
mass density for these runs. Both simulations ran for 
more than seven orbits before merger. As previously 
mentioned, the goal of these simulations is not to provide 
highly accurate models of the mergers but simply to test 
the ID sets: the simulations track the BNS evolution un- 
til the code crashes due to the formation of singularities 



10"' 



10- 



10"' 



A, A , |HC 05 -HC 25I - 


iUia,1 nm 1 6|HC 25 -HC -125I - _ 

mwm 

m n 


i. fl' " 'ft, jn 




- V I 

1 





10 



15 



20 

x/M 



25 30 



35 
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shows the difference between the medium (O.25/M0) and low 
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and no attempt is made to study in detail the forma- 
tion/evolution of the newly created compact object. 

While some of the design parameters of the set BNS80n 
were chosen to produce a meaningful comparison with 
LORENE80 (coordinate separation, total rest mass, zero 
spins), the boost velocity v l was fine-tuned to minimize 
eccentricity. The purpose was to find ID sets less eccen- 
tric than those based on helical symmetry. This idea was 
inspired by the work of Miller [l(| who showed that zero- 
ing the radial velocity could lead to non-negligible orbital 
eccentricities. More recent work done on BBH ID sets 
shows that this eccentricity can be controlled to some ex- 
tent by a careful choice of tangential and transverse veloc- 
ities |56rro0l |. We determined that a boost velocity with 
tangential v$ = 0. 1265c and radial v r = —0.005c com- 
ponents achieves the desired effects for the set BNSSOn. 
The determination of this velocity was done by a correc- 
tive iteration and could be sped up by adopting a more 
efficient technique. Future work will study the applica- 
bility of methods developed for BBH ID sets [H, EH ■ 

Figures [H] and [5] plot the trajectory of one of the NS for 
each simulation and the coordinate separation between 
NS centers respectively. Note that the merger of BNSSOn 
occurs about 0.5 ms earlier than that of LORENE80. 
Figure [TU] shows an estimation of the eccentricity for both 
runs calculated using formula (1) from 61|. Note that by 
the last orbit, the eccentricity exhibited by the BNS80n 
run is about an order of magnitude smaller than that of 
the LORENE80 simulation. This effect can also be ob- 
served in the amplitude of the gravitational waves (GW) 
in Fig. [11] where, for clarity, the signals were shifted in 
time to make the maxima in amplitude coincide. Addi- 
tionally, our ID set seems to possess less "junk" radia- 
tion. Whether this effect persist in the case of BNS with 
generic spins will be studied in future work. 
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FIG. 7: Rest mass density for the Lorene80 (top) and the 
BNS80n (bottom) simulations. 



To further assess the quality of the BNS80n simula- 
tion, we plotted the maximum rest mass density as a 
function of time and compared it with the curve corre- 
sponding to the LORENE80 run (Fig. [12]). Our ID run 
shows oscillations with larger amplitude than those of 
LORENE80 during the first orbits; a behavior already 
discussed in Fig. [T] However, the amplitude diminishes 
with time and, after a couple of orbits (t ~ 13 ms), 
both runs show comparable oscillations. The most im- 
portant question for the assessment of the quality of our 
ID sets is, however, how large are the violations of the 
constraints? Figure [T3] presents a qualitative view of the 
Hamiltonian constraint violation through snapshots that 
compare LORENE80 vs. BNS80n. A quantitative com- 
parison is given by the L2 norm of the Hamiltonian and 
the y component of the momentum constraint violations 
for both runs as a function of time (Fig. [T4l) . As ex- 
pected, at t — the constraint violations in BNSSOn are 
larger than those of LORENE80. However, after a couple 
of orbits, the difference diminishes and the scale of the 
violations for each run become comparable. Other quan- 
titative indicators such as the L\ and L m norms show 
similar behavior. This damping is attributable to the 
well known constraint violation propagation properties of 




FIG. 8: Trajectory of one NS for the runs of Fig. [7] 
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FIG. 9: Coordinate separation between NS centers for the 
runs of Fig. (7) 



the BSSNOK formulation and seems to indicate that this 
violation reduction could possibly be hasten with the use 
of formulations with superior constraint damping char- 
acteristics such as CCZA [H, 63] and ZAc [g4[ (Kastaun 
et al. [34| show CCZ4 simulations where the contraint 
violations fall below BSSNOK levels after only 1ms). A 
study of how these quality control markers depend on the 
grid resolution is important and, due to the correspond- 
ing large computational cost, outside the scope of this 
work. 



B. Spinning binaries evolution 

The most salient characteristic of our method for con- 
structing BNS ID sets is the ability to handle spinning 
stars. To show this, we produced two sets with NS 
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FIG. 10: Eccentricity estimation for the runs of Fig. [7] The 
calculation was done using formula (1) from [611 ]. 
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FIG. 11: Gravitational wave amplitude (bottom) and real 
part (top) of the (2, 2) mode of ^4 for the runs of Fig. [7] Sig- 
nals were shifted to coincide at the amplitude maximum (top 
and bottom) and at the "junk" radiation (lower left inset). 



with spins aligned (BNS80u) and anti-aligned (BNS80d) 
with the orbital angular momentum and evolved them 
from a starting coordinate separation of 80 km through 
their mergers. Both sets have identical NS with rest 
masses Mo = 1.625M Q and dimensionless spin param- 
eters a — J/M 2 ~ 0.33, as detailed in Table HI 

Figure [15] shows the trajectories of one NS for the 
evolution of each set; these curves complement the non- 
spinning BNS (BNS80n) trajectory given in Fig. [8] One 
interesting feature is the presence of the orbital "hang- 
up" [43[. This effect predicts that systems with spins 
aligned with the orbital angular momentum orbit longer 
than those with anti-aligned spins, allowing the shedding 
of excess of angular momentum through the emission of 
gravitational radiation. While the anti-aligned binary 
merged and immediately collapsed to a black hole, the 



FIG. 12: Maximum rest mass density normalized to the initial 
values for the evolution of ID sets LORENE80 and BNS80n 
in Table UJ M represents the ADM mass of the LORENE80 
set. 



aligned case led to the formation of a centrifugally sup- 
ported hypermassive neutron star that survived for about 
20 ms before collapsing. Since the scope of the paper is 
the study of BNS inspirals, no attempts were made to de- 
tect horizons or implement excision-like algorithms that 
would allow for a more accurate modeling of the merger 
and ringdown stages. 

Figures [TTd] and [T7] show the maximum rest mass den- 
sity and the coordinate separation vs. time for the three 
BNS 80 runs. The time to merger of each set depends 
on the values of the boost velocities used for each case 
(reported in Table [l|. We have made an effort to find 
the velocities that would minimize the eccentricity but 
we halted the fine-tuning after reaching what we deemed 
a reasonable precision considering the low resolution of 
our simulations. More accurate evolutions with even 
lower eccentricities are likely to result in merger times 
that are quantitatively different than the values shown 
here. Figure [15] presents the orbital eccentricities of the 
three BNS 80 runs and compares them with that of the 
LORENE80 simulation. The eccentricities achieved dur- 
ing the last orbits of the three BNS80 cases are lower 
than that of the reference LORENE80 simulation. 



V. CONCLUSIONS 

We introduced a new way of constructing initial data 
for binary neutron stars with arbitrary spins and or- 
bital eccentricities. The method only offers approxima- 
tions to the Einstein field equations since, by design, the 
data sets do not satisfy the Hamiltonian and momentum 
constraints. However, by evolving our initial data, we 
showed that after a couple of orbits these constraint vi- 
olations become comparable to those seen in evolutions 
of standard (i.e., irrotational, conformally-flat, helically 
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FIG. 13: Absolute value of the Hamiltonian constraint viola- 
tion for the Lorene80 (top) and the BNSSOn (bottom) sim- 
ulations. The time of the snapshots match those of Fig. [7] 
The color logarithmic ruler is shown at the bottom. 



symmetric, constraint solving) initial data sets. Our 
method consists of a variant of metric superposition that 
addresses two common problems: large stellar shape os- 
cillations and orbital eccentricities. It reduces the former 
to levels similar to those found in current simulations of 
standard sets and offers great control over orbital eccen- 
tricities. Additionally, we see indications that our initial 
data sets possess less "junk" radiation than that found 
in standard sets. 

However, its most important characteristic is the abil- 
ity to handle spinning binaries. We tested our initial 
data by evolving in time initial data for single and bi- 
nary neutron star systems. We showed that our data 
leads to inspirals with orbital eccentricities smaller than 
those seen in standard initial data simulations (up to one 
order of magnitude smaller for the last orbits). We also 
evolved binaries with spins aligned and anti-aligned with 
the orbital angular momentum. The anti-aligned binary 
merges and immediately collapses to a black hole, while 
the aligned case leads to the formation of a centrifugally 
supported hypermassive neutron star that survives for 
several dynamical times before collapsing. 

The work presented here will be followed by studies 
that are outside of the scope of this article due to com- 
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FIG. 15: Trajectory of one NS for the evolution of BNS ID sets 
with spins aligned (BNS80u, left) and anti-aligned (BNS80d, 
right) with the orbital angular momentum. 



putational demands. We plan to explore the viability of 
the "over-correction" technique for binaries with generic 
spins. We will study more efficient ways of selecting the 
stars' boost velocity since the direct trial-and-error ap- 
proach employed here is time and resource consuming. 
We will conduct a more systematical study of the con- 
tent of "junk" radiation in our sets, to find out if the 
reduction in this unwanted quantity is a common feature 
in binaries with arbitrary spins. Another priority is a 
more exhaustive analysis of the characteristics of the ini- 
tial data as a function of coordinate separation. Finally, 
more realistic simulations (higher grid resolutions, dif- 
ferent equations of state, horizon detection and tracking, 
ringdown modeling) and stars with arbitrary spins (which 
would require three-dimensional grids) will be pursued. 



10 



t/M 

500 1000 1500 2000 2500 3000 








5.0 - 






3.0 - 


1.1 




1.0 - 


o" 
II 







1 1.0 

Q. 






Pmax' 







0.9 



20 30 
time [ms] 

FIG. 16: Maximum rest mass density normalized on the initial 
value for simulations of BNS with spins aligned (BNS80u), 
zero (BNS80n) and anti-aligned (BNS80d) with the orbital 
angular momentum. 



t/M 



500 1000 1500 2000 2500 



40 



30 



c\j 20 



10 





























S X. 




\ \ 




\ \ 
\ \ 


BNS80U 


I \ - 
I \ 


BNS80n 


\ V 


BNS80d - 


\ 



CM 

4 § 



10 



20 30 
time [ms] 



40 



FIG. 17: Coordinate separation for the simulations of Fig. 

ED 



This research was supported by the NSF award PHY- 
0855315 and by an allocation of advanced comput- 
ing resources provided by the National Science Foun- 
dation. The computations were performed on Krakcn 
at the National Institute for Computational Sciences 

(http : // www. nics . tennessee.edu / ) . 

t/M 

500 1000 1500 

0.08 



0.06 



c 0.04 

<D 
O 

o 



0.02 



LORENE80 
BNS80U 
BNS80n 
BNS80d 




10 



15 

time [msl 



20 



25 



FIG. 18: Orbital eccentricity for the simulations of Fig. 1161 



3000 4000 



9* 



0.001 



-0.001 



9£ 0.001 
0.000 r*= 



BNS80U 
BNS80n 
BNS80d 




10 20 30 40 50 60 
time [ms] 



Acknowledgments 

It is a pleasure to thank K. Yakunin for useful dis- 
cussions. We are also grateful to S. Bernuzzi, T. Diet- 
rich, F. Galeazzi, T. Font, J. Friedman, W. Tichy and in 
particular L. Rezzolla for comments on the manuscript. 



FIG. 19: Gravitational wave amplitude (bottom) and real 
part (top) of the (2, 2) mode of ^4 for the simulations of Fig. 
ED 



[1] N. Gehrels and P. Meszaros, Science 337, 932 (2012). 
[2] A. Gomboc (2012), 1206.3127. 

[3] B. Sathyaprakash and B. Schutz, Living Rev.Rel. 12, 2 

(2009), 0903.0338. 
[4] J. Abadie et al. (LIGO Scientific Collaboration, Virgo 

Collaboration), Class.Quant.Grav. 27, 173001 (2010), 

astro-ph/1003.2480. 



[5] J. R. Wilson and G. J. Mathews, Phys. Rev. Lett. 75, 
4161 (1995). 

[6] J. R. Wilson, G. J. Mathews, and P. Marronetti, Phys. 

Rev. D54, 1317 (1996), gr-qc/9601017. 
[7] P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963). 
[8] S. Oslowski, T. Bulik, D. Gondek-Rosinska, and K. Bel- 

czynski (2009), 0903.3538. 



11 



I. Kowalska, T. Bulik, K. Belczynski, M. Dominik, 
and D. Gondek-Rosinska, Astron.Astrophys. 527, A70 
(2011), 1010.0511. 

M. Miller, Phys. Rev. D 69, 124013 (2004), gr- 
qc/0305024. 

T. W. Baumgarte, G. B. Cook, M. A. Scheel, S. L. 
Shapiro, and S. A. Teukolsky, Phys. Rev. Lett. 79, 1182 

(1997) , gr-qc/9704024. 

T. W. Baumgarte, G. B. Cook, M. A. Scheel, S. L. 
Shapiro, and S. A. Teukolsky, Phys. Rev. D 57, 6181 

(1998) . 

L. Bildsten and C. Cutler, Astrophys. J. 400, 175 (1992). 
C. Kochanek, Astrophys. J. 398, 234 (1992). 
S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. 
Rev. D 56, 7740 (1997). 

H. Asada, Phys. Rev. D 57, 7292 (1998), U RL 

|http : //link . aps . org/doi/ 10 . 1 103/ Phy sRevD . 57 . 7292 

S. Teukolsky, Astrophys. J. 504, 442 (1998). 

M. Shibata, Phys. Rev. D 58, 024012 (1998). 

P. Marronetti, G. J. Mathews, and J. R. Wilson, Phys. 

Rev. D58, 107503 (1998), gr-qc/9803093. 

S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys. 

Rev. Lett. 82, 892 (1999), gr-qc/9810072. 

P. Marronetti, G. J. Mathews, and J. R. Wilson, Phys. 

Rev. D 60, 087301 (1999), arXiv:gr-qc/9906088. 

K. Uryu and Y. Eriguchi, Phys. Rev. D61, 124023 

(2000) , gr-qc/9908059. 

K. Uryu, M. Shibata, and Y. Eriguchi, Phys. Rev. D62, 

104015 (2000). 

E. Gourgoulhon, P. Grandclement, K. Taniguchi, 
J. Marck, and S. Bonazzola, Phys. Rev. D 63, 064029 

(2001) . 

K. Taniguchi, E. Gourgoulhon, and S. Bonazzola, Phys. 
Rev. D64, 064012 (2001), gr-qc/0103041. 
K. Taniguchi and E. Gourgoulhon, Phys. Rev. D66, 
104019 (2002), gr-qc/0207098. 

J. A. Faber and F. A. Rasio, Living Rev.Rel. 15, 8 (2012), 
1204.3858. 

K. Uryu, F. Limousin, J. L. Friedman, E. Gourgoulhon, 
and M. Shibata, Phys.Rev.Lett. 97, 171101 (2006), gr- 
qc/0511136. 

G. B. Cook and T. W. Baumgarte, Phys. Rev. D 78, 

104016 (2008). 

K. Uryu, F. Limousin, J. L. Friedman, E. Gourgoul- 
hon, and M. Shibata, Phys.Rev. D80, 124004 (2009), 
0908.0579. 

M. Anderson, E. W. Hirschmann, L. Lehner, S. L. 
Liebling, P. M. Motl, et al., Phys.Rev. D77, 024006 
(2008), 0708.2720. 

R. Gold, S. Bernuzzi, M. Thierfelder, B. Brugmann, and 

F. Pretorius, Phys.Rev. D86, 121501 (2012), 1109.5128. 
W. E. East, F. M. Ramazanoglu, and F. Pretorius, 
Phys.Rev. D86, 104053 (2012), 1208.3473. 

W. Kastaun, F. Galeazzi, D. Alic, L. Rezzolla, and J. A. 
Font (2013), 1301.7348. 

A. Lyne, M. Burgay, M. Kramer, A. Possenti, R. Manch- 
ester, et al., Science 303, 1153 (2004), astro-ph/0401086. 
W. Tichy, Phys.Rev. D84 024041 (2011), URL 
|10.1103 /PhysRevD .84.024041] 

P. Marronetti and S. L. Shapiro, Phys. Rev. D68, 104024 
(2003), gr-qc/0306075. 

T. W. Baumgarte and S. L. Shapiro, Phys.Rev. D80, 

064009 (2009), 0909.0952. 

W. Tichy, Phys.Rev. D86, 064024 (2012). 



[40] URL http: //www . lorene . obspm. f r/ 
[41] URL http : / /einsteintoolkit . org/ 
[42] F. Loffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, 

R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnet- 

ter, et al, Class. Quant. Grav. 29, 115001 (2012). 
[43] M. Campanelli, C. Lousto, and Y. Zlochower, Phys.Rev. 

D74, 041501 (2006), gr-qc/0604012. 
[44] N. Stergioulas, Living Reviews in Relativity 6 (2003). 
[45] N. Stergioulas and J. L. Friedman, Astrophys. J. 444, 

306 ( 1995). 

[46] URL http://www.gravity.phys.uwm.edu/rns/ 

[47] H. Komatsu, Y. Eriguchi, and I. Hachisu, Mon. Not. R. 

Astron. Soc. 237, 355 (1989). 
[48] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro- 
phys. J. 422, 227 (1994). 
[49] R. A. Matzner, M. F. Huq, and D. Shoemaker, Phys. 

Rev. D 59, 024015 (1999). 
[50] P. Marronetti, M. F. Huq, P. Laguna, L. Lehner, R. A. 

Matzner, and D. Shoemaker, Phys. Rev. D 62, 024017 

(2000), gr-qc/0001077. 
[51] E. Schnetter, F. Herrmann, and D. Pollney, Phys. Rev. 

D 71, 044033 (2005), gr-qc/0410081. 
[52] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. 

Phys. Suppl. 90, 1 (1987). 
[53] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, 

024007 (1998), gr-qc/9810065. 
[54] J. D. Brown, P. Diener, O. Sarbach, E. Schnetter, and 

M. Tiglio, Phys.Rev. D79, 044023 (2009), 0809.3533. 
[55] L. Baiotti, I. Hawke, P. J. Montero, F. Loffler, L. Rez- 
zolla, N. Stergioulas, J. A. Font, and E. Seidel, Phys. 

Rev. D 71, 024035 (2005), gr-qc/0403029. 
[56] S. Husa, M. Hannam, J. A. Gonzalez, U. Sperhake, 

and B. Bruegmann, Phys.Rev. D77, 044037 (2008), 

0706.0904. 

[57] B. Walther, B. Bruegmann, and D. Mueller, Phys.Rev. 

D79, 124040 (2009), 0901.0993. 
[58] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom, 

G. Lovelace, et al., Class.Quant.Grav. 24, S59 (2007), 

gr-qc/0702106. 

[59] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. 

Pfeiffer, et al., Phys.Rev. D76, 124038 (2007), 0710.0158. 
[60] A. H. Mroue, H. P. Pfeiffer, L. E. Kidder, and S. A. 

Teukolsky, Phys.Rev. D82, 124016 (2010), 1004.4697. 
[61] W. Tichy and P. Marronetti, Phys.Rev. D83, 024012 

(2011), 1010.2936. 
[62] S. Bernuzzi and D. Hilditch, Phys.Rev. D81, 084003 

(2010), 0912.2920. 
[63] D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, 

and C. Palenzuela, Phys.Rev. D85, 064040 (2012), 

1106.2254. 

[64] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, 

W. Tichy, et al. (2012), 1212.2901. 
[65] Other constraints could also be present if additional fields 

(such as electromagnetic) are included. 
[66] During the preparation of this article, Kastaun et al. [H|] 

confirmed that this effect is even more pronounced when 

evolving the sets with the CCZ4 formulation instead of 

BSSNOK. 

[67] This is achieved by the use of a conformal factor as an 
auxiliary field derived from solving the Hamiltonian con- 
straint, and by corrections to the superposed shift vector 
that solve the momentum constraint [33J. 



